#Ellipse and the Malinvaud example

pellipse <- function(b, c, df, level, xlab = "", ylab = "",add=FALSE)

{

        d <- sqrt(diag(c))

        df <- c(2, df)

        phase <- acos(c[1, 2]/(d[1] * d[2]))

        angles <- seq(-(pi), pi, len = 51)

        mult <- sqrt(df[1] * qf(level, df[1], df[2]))

        xpts <- b[1] + d[1] * mult * cos(angles)

        ypts <- b[2] + d[2] * mult * cos(angles + phase)

	if(add) lines(xpts, ypts)

	else plot(xpts, ypts, type = "l", xlab = xlab, ylab = ylab)

}



b1 <- c(0.133, 0.55)



b2 <- c(-0.021, 0.235)



C1 <- matrix(c(3.9e-05, -0.000147, -0.000147, 0.01219) , nrow = 2, ncol = 2)



C2 <- matrix(c(0.0025, -0.0038, -0.0038, 0.0059) , nrow = 2, ncol = 2)



postscript("malinvaud.ps",width=7, height=5, horizontal=FALSE)

par(mfrow=c(1,2))

pellipse(b1,C1,14,.95,xlab="b.gdp",ylab="b.investment")

title("95% Confidence Ellipse","GDP and Investment on Imports")

pellipse(b2,C2,14,.95,xlab="b.gdp",ylab="b.consumption")

title("95% Confidence Ellipse",  "GDP and Consumption on Imports")

dev.off()